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ABSTRACT 

A  Green's  function  method  is  used  to  find  the  low  temperature 
change  in  the  specific  heat  due  to  a  (110)  surface  on  a  simple 
cubic  monatomic  lattice.  Two  separate  first  neighbor  force  con¬ 
stant  models  are  used  for  the  calculation:  the  first  assumes  that 
the  atomic  motion  normal  to  the  surface  is  uncoupled  from  motion 
parallel  to  the  surface;  the  second  is  the  familiar  two  force  con¬ 
stant  model  popularized  by  Montroll  and  Potts  .  Both  models  are 
anisotropic  in  the  surface  and  neither  satisfies  the  condition  of 
rotational  invariance.  Analytic  expressions  are  found  for  the  sur¬ 
face  mode  dispersion  relations  and  for  the  low  temperature  specific 
heat.  It  is  found  that  for  small  deviations  from  isotropy,  the 
change  in  the  specific  heat  is  independent  of  the  model  and  is  the 
same  for  the  (110)  and  (100)  surfaces. 
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13  ABSTRACT 


A  Green's  function  method  is  used  to  find  tne  low 
temperature  cnange  in  the  specific  neat  due  to  a  (110) 
surface  on  a  simple  cubic  monatomic  lattice.  Two  separate 
first  neignbor  force  constant  models  are  used  for  the  cal¬ 
culation:  tne  first  assumes  that  tne  atomic  motion  normal 
to  tne  surface  is  uncoupled  from  motion  parallel  to  tne 
surface;  tne  second  is  the  familiar  two  force  constant  model 
popularized  by  Montroll  and  Potts.  Both  models  are  aniso¬ 
tropic  in  tne  surface  and  neither  satisfies  the  condition 
of  rotational  invariance.  Analytic  expressions  are  found 
for  the  surface  mode  dispersion  relations  and  for  the  low 
temperature  specific  neat.  It  is  found  that  for  small 
deviations  from  isotropy,  the  change  in  the  specific  neat 
is  independent  of  the  model  and  is  the  same  for  the  (110) 
and  (100)  surfaces. 


DD,FrJ473 

S/N  0 10 )•  807. 6  80 t 


( P  A'iL  1] 


Unclassif ied 


Si  L  dfi  t\  ( 


2 


THE  SPECIFIC  HEAT  OF  AN  ANISOTROPIC  SURFACE 

There  have  been  a  number  of  studies  of  the  change  in  the 
specific  heat  due  to  the  presence  of  a  surface  from  both  the 
continuum  and  the  discrete  lattice  points  of  view.'1'  In  every 
study  for  which  an  analytic  result  is  obtained,  however,  the  cal¬ 
culation  has  been  restricted  to  force  models  that  are  either  iso¬ 
tropic  in  the  balk  of  the  crystal  or  isotropic  in  the  plane  of  the 
surface. 

In  this  paper,  we  will  obtain  the  change  in  the  specific  heat 
due  to  a  (110)  surface  bn  a  simple  cubic  monatomic  lattice.  The 
lattice  will  be  described  by  two  models,  both  of  which  use  two  first 
neighbor  force  constants.  The  major  interest  in  these  calculations 
comes  from  the  fact  that  the  two  models  used  are  anisotropic  in  the 
(110)  surface  and  from  the  fact  that  the  models  are  sufficiently 
simple  for  the  entire  calculation  to  be  carried  out  analytically. 

The  models  have  the  disadvantage,  however,  that  they  do  not  satisfy 
rotational  invariance  and  therefore  the  results  must  be  interpreted 
accordingly . 
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I.  SURFACE  MODES 

.  The  Cartesian  coordinate  system  will  be  situated  s.  .  - Va  that  the 
x  and  y  unit  vectors  lie  in  the  plane  of  the  surface  (?* a  Figure  1) . 
In  the  harmonic  approximation,  the  force  on  the  (r  n  atv-M  is: 

V0)  =  mUa(0)  ^(O^u  U)  (1.1) 

where  m  is  the  mass  of  each  atom  and  u  (/.)  is  the  dsplacement 

ot 

i.  L. 

from  equilibrium  of  the  t  atom  in  the  cr-Cartesian  direction. 

The  particular  crystal  models  to  be  considered  are  determined  by 
specifying  the  values  of  the  force  constants  §^^(01).  In  this 
paper,  the  force  constants  <f>ag((H)  are  zero  unless  the  index  l 
refers  to  one  of  the  six  nearest  neighbors  of  the  0l  atom.  The 
labelling  we  will  use  for  four  of  these  neighbors  is  shown  in 
Figure  1.  Atoms  labelled  5  and  6  are  located  at  ±a  (1,0,0). 

A.  Model  I ,  Equations  of  Motion 

The  first  model  to  be  considered  involves  central  and 
non-central  nearest  neighbor  forces  defined  such  that  the  motion 
of  the  ions  in  any  one  of  the  Cartesian  directions  is  uncoupled 
from  the  motion  in  any  other  direction.  The  equations  of  motion 
for  a  bulk  atom  in  this  irodel  are: 

mU  (  0)  =  6[u  ( 5)  +u  (6)  -2u  (  0)  ]  +(3[u  ( 1)  +u  ( 2)  +u  (3)  +uv ( 4)  -4uv (  0)  ] 
mUy  (0 )  =  p  [  uy  ( 5)  +uy  ( 6)  -2uy  (0  )  ]  +or[uy  ( 1)  +uy  (  2)  +uy  ( 3)  +uy  ( 4)  -4uy  { 0 )  ] 
mUzt0)  “  3[u2(5)+uz(6)-2uz(0)]+a[uz(l)+uz(2)+uz(3)+uz(4)-4uz(0)] 


(1.2) 


The  coefficient  g  is  the  non-central  force  constant  defining  the 
restoring  force  when  the  atoms  in  the  bond  are  displaced  perpen¬ 
dicular  to  the  bond  direction.  The  coefficient  6  is  the  central 
force  constant  for  the  bonds  extending  in  the  x-direction.  The 
force  constants  6  and  g  are  related  to  the  elastic  constants  by: 

6  *  aoCll 

(1.3) 

9  "  ao°44=  -  aDC12 

and  the  force  constant  a  is: 

a  -  \  (6  +  g)  (1.4) 

These  relations  result  from  identifying  the  equations  of  motion  in 

the  long  wavelength  limit  with  the  Christoffel  equations  of  elasticity. 

The  interest  in  this  very  simple  model  stems  from  the  fact 

that  it  giver?  rise  co  Rayleigh  waves  that  are  split  off  from  the 

bulk  modes  even  though  the  motion  of  the  ions  normal  to  the  surface 

is  completely  uncoupled  from  the  motion  of  the  ions  parallel  to  the 

2  3 

surface.  This  supports  the  previous  observation  ’  that  a  suffi¬ 
cient  criterion  for  the  presence  of  Rayleigh  waves  is  that  the  cut 
bonds  be  oblique  to  the  sarface. 

B.  Model  II,  Equations  of  Motion 

The  second  model  is  the  now  familiar  model  first  suggested 

4  5 

by  Rosenstock  and  Newell  (and  made  popular  by  Montroll  and  Potts  ) . 

This  model  uses  two  force  constants:  a  central  force  constant  y 
which  bives  ‘the  restoring  force  when  the  atoms  in  a  bond  are  dis¬ 
placed  in  the  direction  of  the  bond,  and  a  non-central  force  constant 
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X  which  gives  the  restoring*  force  when  the  displacements  are  per¬ 
pendicular  to  the  bond.  Due  to  our  choice  of  coordinate  axes, 
these  force  constants  lead  to  a  coupling  of  the  equations  of 
motion  in  the  y  and  z  direction.  The  equations  of  motion  for  a 
bulk  atom  in  this  model  are: 

mlix(0)  =  y[ux(5)+ux(6) -2ux(0)  ]  +  X[ux(l)+ux(2)+ux(3)+ux(4)  -4ux(0)  ] 
mUy(0)  «  X[uy(5)+uy(6)-2uy(0)]  +  |(y-\)[uz(1)-uz(2)+uz(3)-uz(4)] 

+  i(Y+X)[uy<l)+uy(2)+uy(3)+uy(4)-4uy(0)] 
mUz(0)  =  X[uz(5)+uz(6) -2uz(0)  ]+  |(y“X)  [uy ( 1)  -uy (2)  +uy (3)  -uy (4)  ] 

+  £(y+X)  [uz(l)+uz(2)+uz(3)+uz(4)  -4uz(0)  ]  (1.5) 

The  force  constants  and  elastic  constants  are  related  by: 

Y  =  aoCll 

(1.6) 

X  ”,  aoC44  ~  "  aoC12 

The  interest  in  this  model  is  partly  pedagogical  and  partly 

due  to  the  ability  to  compare  our  specific  heat  results  with  the 

X  6 

same  model  calculation  for  the  (100)  surface.  ’ 

We  assume  for  both  models  that  the  displacements  of  the  atoms 
can  be  written  as: 

u  (t)  =  U  exp[ik*x(t)+  iuut]  (1.7) 

where  uu  is  the  angular  frequency  and  k  is  the  wave  vector.  To 
obtain  the  dispersion  relations  for  the  bulk  modes,  we  make  the 
usual  assumption  of  cyclic  boundary  conditions  in  the  three 
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Cartesian  directions. 

By  substituting  Equation  (1.?;  into  the  equations  of  motion 
for  the  two  models,  we  obtain  matrix  equations  of  the  form 

(Model  I) 


D(I>-m.2 


Uy  -  0 


(1.8) 


(Model  II) 


D(II)-  muu2 

yy 


D(II>-m.2 

yy 


(1.9) 


It  is  easily  seen  that  the  eigenvectors  which  diagonalize  the 
above  matrices  are  independent  of  the  value  of  the  wave  vector 
and  are  given  by: 


(Model  I) 


e^(kj)  =  0 


(Model  II) 
„(II) 


(kj)  - 


(1.10) 


1  •  -1 


(1.11) 


r 


I- 


*fi5TO"»y*^ 


The  diagonalized  matrix  elements  when  set  equal  to  zero  give  the 
frequencies  of  the  bulk  modes.  These  are: 


(Model  I) 


m-rj(k)  =  =  SD^^+pDgCk) 


ma)^3(k)  =  Dyy}  =  PD^^+o-DgCk) 


(1.12) 


(Model  II) 


where : 


mu)^(k)  =  D^I}  =  YD1(k)+XD2(k) 

•4v  -  °r  -  -C 

'«  XD1(k)+|(Y+X)D2(k)-|(Y-X)D3(k)  (1.13) 

-ft)  -  D<»>  +  d”*> 

«  XD1(k)+|(Y+X)D2(k)+|(Y-X)D3(k) 


D,(k)  =  2(1  -  cos  k  a  ) 

JL  ~  X  O 

D2(k)  =  4(1  -  cos  2  kya-^cos  \  k^) 
Dg(k)  =  4  sin  \  sin  \  ^zai 


and 


a,  =  \T2f  a 
x  o 


(1.14) 

(1.15) 


For  comparison  with  the  surface  modes,  we  set  the  component 
kz=  0  to  obtain  the  dispersion  relations  for  the  bulk  modes  that 
propagate  parallel  to  the  surface.  This  gives: 

(Model  I ,  bulk) 

j? C k  k  )  =  26(1  -  cos  k  a  )+  4j3(]  -  cos  \  k  a,) 

1  x  y  x  o  K  z  y  1 


muj. 


(1.16) 


ra'jjg  3(kxky)  =  23(1  -  cos  kxao)+  4a(l  -  cos  \  kyai) 


I 
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(Model  II,  bulk) 

A 

mu^(kxky)  =  2y(l  -  cos  kxaQ)  +  4X(1  -  cos  |  k^) 

muj?  0(k  k  )  =  2X(1  -  cos  k  a  )  +  2(y+\)(l  -  cos  l  k  a .) 

*  a  y  aw  y 

(1.17) 

These  dispersion  relations  are  shown  in  Figure  2. 

The  models  we  have  chosen  are  sufficiently  simple  that  the 
dispersion  relations  for  the  surface  modes  can  be  determined  by 

n 

either  a  direct  boundary  value  calculation  or  a  Green's  function 

g 

approach.  The  mathematical  steps  for  the  Green's  function  method 
are  the  more  cumbersome,  however,  and  therefore  we  will  defer  to 
the  boundary  value  method  in  this  section. 

As  indicated  in  Figure  1,  we  will  assume  the  crystal  occupies 
the  half-space  where  z  >  0.  Thus,  we  will  assume  periodic  boundary 
conditions  in  the  x  and  y  directions,  but  replace  the  boundary  con¬ 
dition  in  the  z  direction  by  the  requirement  that  the  sum  of  all 
forces  between  the  atoms  on  opposite  sides  of  the  surface  must  van¬ 
ish.  For  our  models,  these  are  the  forces  on  atom  0  due  to  atoms 
1  and  2, 

C.  Model  I,  Boundary  Conditions 

The  boundary  conditions  for  Model  x  can  be  written: 

Fx  =  °  =  g[ux(lUux(2)-2ux(0)] 

Fy  =  0  «  olu  (l)+uy(2)-2uy(0)] 

Fz  -  0  =  a[u2,(l)+uz(2)-2u;?(0)] 


(1.18) 


Since  the  x,  y,  and  z  atomic  displacements  in  this  model  are  inde¬ 
pendent,  these  boundary  conditions  are  satisfied  by  displacements 
of  the  form  (1.7).  Upon  substitution,  we  obtain  the  same  complex 


equation  for  each  condition  in  (1.18);  namely: 

i|k  a, 

e  =  cos  2  kyai 


(1.19) 


By  writing  k  in  the  form: 

Z 

kz  =  kz  +  ikz  (1.20) 

reduces  to: 

(1.21) 

in  Figure  3.  The  imagin¬ 
ary  component  k^  (called  the  attenuation  coefficient)  leads  to  an 
amplitude  of  vibration  in  Equation  (1.7)  which  varies  exponen¬ 
tially  with  increasing  distance  from  the  surface.  Solutions  wheve 
k*  is  positive  are  classified  as  surface  waves. 

Z 

The  dispersion  curves  for  the  surface  modes  are  found  by  sub¬ 
stituting  the  value  of  the  wave  vector  component  k  imposed  by  the 

z 

boundary  conditions  (1.18)  into  Equations  (1.12).  This  gives  yhe 
simple  result: 

(Model  I,  surface) 
o 

muu,  (k  k  )  =  26(1  -  cos  k  a  )+  3(1  -  cos  k  a.) 

-*•  x  y  x  o  y  j. 

2 

nut) q  o(k  k  )  =  23(1  -  cos  k  a  )  +  a(l  -  cos  k  a,) 
z j ^  y  xo  y  i 


where  k  and  k  are  real,  Equation  (1.19) 

Z  Z 


kz  ai  =  0 


"kz  al 


=  2(1  +  COS  k  ap 


These  relationships  are  shown  graphically 


(1.22) 
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Thefc*’  dispersion  curves  are  shown  in  Figure  2.  Of  interest  is 

to  note  that  for  small  values  of  the  wave  vector  k,  the  dispersion 

curves  for  the  surface  modes  and  for  the  bulk  modes  are  the  same 
4 

to  order  k  .  Hen^e,  the  velocities  of  sound  for  the  surface  modes 
are  the  same  as  the  velocities  of  scu-  d  for  the  bulk  modes  having 
the  same  polarization  and  propagation  direction. 

D.  Model  II,  Boundary  Conditions 

The  boundary  conditions  for  Model  II  can  be  written: 

Fx  =  0  =  X[ux(l)+ux(2)-2ux(0)] 

Fy  -  0  =  £(Y+X)[uy(l)+uy(2)-2uy(0>]+  £(y-X)  [uz(l) -uz(2)  ] 

Fz  =  0  =  £(y+X)  [uz(1)+uz(2)  -2uz(0)]+  i(Y-X)  [uy (1) -uy (2)  ] 

(1.23) 

For  this  model,  x-displacements  of  the  atoms  are  independent  of 
the  y  and  z  displacements.  By  substituting  displacements  of  the 
form  (1.7)  into  the  equation  for  the  x  component  of  the  force,  we 
obtain  Equation  (1.19)  relating  the  wave  vector  components  ky  and 
kz.  For  this  particular  mode,  the  steps  are  exactly  the  same  as 
those  for  Model  I,  and  hence  we  obtain  for  the  surface  mode: 

(Model  II,  surface) 

rau^  (k  k  )  =  2y(1  -  cos  k  a  )+  X(1  -  cos  k  a,)  (1.24) 

x  x  y  x  o  y  x 

[Jpon  substituting  Equation  (1.7)  into  the  boundary  conditions 
for  the  y  and  z  component  of  the  force,  we  obtain  two  coupled  equa¬ 
tions  in  the  Fourier  coefficients  U  and  U„.  In  matrix  form  these 

y  z 

are: 
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and 

tan  \  =  ±  tan  *  kyal  (1.28) 

In  expressions  (1.26)  and  (1.28),  the  +  (-)  sign  belongs  with  the 
mode  j  =  2  (3) . 

Ordinarily  when  two  components  of  the  motion  are  coupled, 
the  form  of  the  displacements  which  must  be  used  to  satisfy  the 
boundary  conditions  is  not  Equation  (1.7)  but  rather  a  linear  com¬ 
bination  of  the  type: 

u  ( l )  =  £  exp[ik^  •x(-t)+iwt]  (1.29) 

a  •  n  o  Of  ~ 

J  a  y  O 

where  the  label  j  is  added  since  the  quantities  U  and  k  may  differ 

Of  ~ 

for  the  two  modes.  This  procedure  has  been  fully  illustrated  by 
Gazis,  Herman,  and  Wallis.  The  models  used  in  this  paper,  towever, 
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are  sufficiently  simple  that  the  equation  that  results  from  this 
linear  combination  is  a  perfect  square  of  the  determinant  of  the 
matrix  in  Equation  (1.25) .  Thus  the  straightforward  calculation 
shown  above  gives  the  same  results  for  the  two  coupled  modes  as 
the  more  involved  method. 

Since  the  real  part  of  kz  is  not  zero  for  these  two  modes, 

they  are  properly  classified  as  generalized  Rayleigh  waves.  The 

dependence  of  k  upon  k  is  shown  in  Figure  3.  Of  interest  is 
^  y 

the  fact  that  the  attenuation  coefficient  k1  increases  with  k  to 

z  y 

a  finite  value  at  the  zone  boundary  whereas  for  the  previous  cases, 
the  attenuation  coefficient  increased  without  bound.  For  the 

vibration?...  modes  that  have  an  infinite  attentuation  coefficient, 
only  the  atom,  in  the  surface  layer  are  displaced  from  equilibrium. 

Tne  surface  mode  frequencies  are  obtained  by  substituting  the 
condition  (1,26)  into  the  frequency  expressions  (1,13)  for  the  bulk 

modes.  The  result  is  two  degenerate  modes  whose  frequencies  are 
given  by: 

(Model  II,  surface) 

"'",2,3<bxky)  ”  2*O-C0S  kxa0)4  (1-cos  kyap  O  .30) 

The  polarizations  of  the  degenerate  modes  for  both  models  are  given 
by  the  eigenvector  components  in  (1.10)  and  (1.11).  An  unusual 

feature  of  both  models  is  the  fact  that  we  obtain  three  separate 
surface  mode  branches  corresponding  to  the  three  diffei’ent  polar¬ 
izations  the  atomic  motion.  In  continuum  calculations,  for 
instance,  only  one  Rayleigh  wave  is  obtained. 


13 


The  surface  mode  dispersion  relations  for  Model  II  have  been 
9 

found  previously  and  are  shown  in  Figure  2.  The  solution  for 
Model  I  and  Model  II  are  identical  in  the  isotropic  limit  when 
Y=X  and  ^=j3,  tut  they  diverge  for  other  values  of  the  force  constants. 


II.  THE  SPECIFIC  HEAT 

As  in  the  surface  mode  calculation,  the  low  temperature  spec¬ 
ific  heat  due  to  the  surface  perturbation  can  be  found  in  two  ways. 
One  method  deals  with  the  explicit  determination  of  the  surface 

inodes  and  the  perturbed  bulk  modes  prior  to  calculating  the  specific 
3 

heat.  The  second  method  sidesteps  the  problem  of  determining  the 
perturbed  modes  by  using  a  Green's  function  approach.  Since  the 
Green's  function  method  is  particularly  convenient  for  low  temper¬ 
ature  calculations  v/here  the  long  wavelength  expressions  for  fre¬ 
quency  can  be  used,  we  will  apply  this  second  method  to  our  present 
problem. 

The  method  we  will  follow  is  that  reported  by  Maradudin  and 

Ashkin”^  and  presented  fully  by  Maradudin  and  Wallis ^  (hereafter 

o 

referred  to  ’  MW) .  The  method  uses  the  function  n(y  )  which  is 
defined  as 


n(y2)  =  S  ~2~~2 — 

kj  (  y2+w^(k) 


y2+uu2(k)  j 


(2.1) 


where  w.(k)  and  uu.(k)  are  the  frequencies  of  the  normal  modes  of 
J  ~  J  ~ 

the  crystal  with  a  surface  and  without  a  surface,  respectively. 

2 

MW  have  shown  that  when  the  function  fi(y  )  has  a  logarithmic  sing 
ularity  in  the  limit  jy|-  0, 


lim 

iy|~  0 


0(y2) 


~A  logjyl  +  o(  log  j  y  | ) 


(2.2) 
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then  the  change  in  the  specific  heat,  £Cy,  is  given  by: 
ACy(T)  =  6A  £(3)kB(kBT/ft)2  +  o(T2) 


(2.3) 


where  CO)  is  the  Riemann  zeta  function,  kg  is  Boltzraan's  constant, 
and  T  is  the  absolute  temperature.  The  specific  heat  problem  then 
reduces  to  the  problem  of  determining  the  coefficient  A  in  (2.2). 

MW  have  shown  the  geneval  form  for  the  function  Q(y2)  to  be: 


fi(y2)  -  -  S 

kj 


tOy  >-kj,-y2) 

Cy2+^(k)]2 

j  ~ 


(2.4) 


where  the  term  t(kj  , -kj  , -y2)  is  the  solution  of  the  integral 
equation : 


tOj,  k'j’,-y2)  =  V(kj,  k*  j  ‘  ,-y2) 


k"j”  y2+u)j„(k") 


t(k"j" 


(2.5) 


and  where: 


V(kj  ,  -k  'j  ' ,  -y2)  =  ji-  6(k>kj6(k,',-kj 


x  x 


y  y 


e“<~J>,|’cs(0Wes(‘£'j')6tz+i 

x  [1  -  e^'S^Hl  -  e_i~ # '^(^ 3 


(2.6) 


Here,  |  L,^  is  the  distance  between  the  two  crystal  surfaces, 


The  discussion  in  this  section,  then,  will  proceed  to: 

2 

(1)  evaluate  the  term  V(kj , -k 1 j ’ , -y  )  for  the  two  models,  (2) 

2 

solve  the  integral  equation  in  (2.5)  for  t(kj,k’j’,-y  ),  (3) 

/v  /w 

perform  the  integration  over  the  Brillouin  zone  in  the  long 
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wavelength  limit  indicated  in  Eq.(2.4),  and  (4)  obtain  an  analytic 
expression  for  the  coefficient  A  for  both  models  by  taking  the 
limit  in  Eq. (2. 2) . 

If  the  eigenvectors  e^(kj)  are  independent  of  the  wave  vector 
k,  then  for  a  particular  value  of  the  index  the  sum  over  a  and 
P  in  Eq. (2. 6)  can  be  carried  out  give: 

S  ea(kj)(|>ap(Ot)eg(-k/j/)  -  E  ea(j)(i>a3(0t)ep(j/) 

Of  p  ^  P 

-(2-7> 

This  result  is  easily  seen  to  be  true  since  the  ^^(0-6)  matrix  has 
exactly  the  same  form  as  the  dynamical  matrix  for  a  given  choice 
of  l. 

The  number  of  values  of  the  index  -t  is  determined  by  the 
number  of  bonds  in  the  perfect  crystal  that  cross  the  boundary. 

For  our  models,  this  number  is  two  representing  the  bonds  between 
atoms  0—1  and  0-2.  Thus  the  terras  Cj(t)  become: 

(Model  I) 

C-^1)  =  Cx(2)  =  p 

C2(1)  =  C2(2)  =  a 

C3(l)  =  C3(2)  =  a  (2-8) 

(Model  II) 

C^l)  =  Cx(2)  =  X 

C2(i)  =  X  ;  c2(2)  =  y 

C3(l)  =  y  ;  C3(2)  =  X  (2.9) 


Equation  (2.6)  can  now  be  written: 


V(kj,-k'.i',-y2)  -  6j,i'6<kx-|1x>6<kPky)/flj2v-t(~J)v'<'~''j')  <2- 
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where: 

v.tQ5j) 


=  2 


C.U ) 

J-L 


mL 


z 


eiik-xU) 


sin|k' 


x(t) 


(2.11) 


Being  able  to  express  the  V(kj , -k '  j  ' , -y  )  term  in  separable 
form  as  in  (2.10)  is  an  important  step  in  the  evaluation  of  the 
.specific  heat.  The  ability  to  do  this  led  MW  to  use  a  model  which 
contained  only  central  interactions.  In  our  calculation,  the  fact 
that  the  eigenvectors  are  independent  of  the  wave  vector  k  allows 


us  to  write  (2.10)  for  a  wide  range  of  different  models. 

By  substituting  the  expression  (2.10)  into  Equation  (2.5), 
we  obtain  the  equation: 


t(kj,k’j\-y2) 


-  £ 
l 


£ 

k"j 


6 (k  -k") 6 (k  -k") 
x  x'  v  y  y ' 


■  (“k"j") 

y2+u;|"(k") 


t(k"j" 


(2.12) 


The  solution  of  this  integral  equation  with  separable  kerne]  is 
straightforward  with  the  result: 


-1 


t(kj,-kj,-y  )  =£/V^(kj)CI-M(kxkyj,-y2)]u,v^(-kj) 


It 


(2.13) 


where: 


2  ^  v » (~kj) v  /(kj) 

"tt'(kxky^-y >  -  -  s  -a  "a*  — 

k  y  +^(k) 

Z  J  ~ 


(2.14) 


Finally,  the  function  Q(y)  in  Equation  (2.4)  becomes: 


Q(y2)  -  -  £  £ 
kj  If 


v/(kj)v^,(-kj) 
[y2+ju2(k)  ]2 


[ 


I-M(kxkyj,-y2)] 


-1 


a 1 


(2.15) 
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Strictly  speaking,  our  choice  of  a  non-primitive  unit  cell 

for  this  problem  leads  to  six  phonon  branches  it  the  first  Brilloui 

zone  as  shown  in  Figure  2.  The  optical  phonon  branches,  however, 

are  extensions  of  the  acoustic  phonon  branches  from  the  neighboring 

zone.  This  means  that  by  using  an  extended  zone  scheme  we  can 

treat  the  problem  as  if  there  are  only  three  independent  phonon 

branches.  Therefore  the  sum  over  "wave  vector  k  and  branch  index 

j  can  be  converted  to  integrals  as: 

ir/a  TT/a,  2rr/a, 

g  v  i.  x  3 

££  -  SP_z-!l  /  dk  f  dk  /  dk  S  (2.16) 

k  j=l  32tt  J  X  J  y  J  j=l 

-rr/ao  -rr/a-^  -2^/^ 

where  SQ  is  the  surface  area  of  the  crystal  after  the  cut  is  made, 

and  the  limits  on  the  integral  over  k  indicate  the  use  of  the 

z 

extended  zone  scheme. 

tVe  separate  out  the  sum  over  j  from  Equation  (2.15)  by 
defining:  2 

n(y2)  -S  Q.(y2)  (2.17) 

j“l  3 


In  this  expression  we  have  obsorbed  the  integral  over  k  into  the 

z 

o 

term  J^'(kxkyj>~y  )  as: 


JU'(kxV>‘y2) 


v, (kj) v^  /(-kj) 
Cy2+oj2(k)  ] 2 

J 


(2.39) 
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In  terms  of  integrals  the  M(k  k  J,-y  )  matrix  becomes: 

x  y 


2rr/a , 


_  La 

M, .  /(k  k  , -y  )  =  -  —r 
ll  v  x  y’  J  7  4rr 


Vi  r 

4rr  J 


dk. 


z  [y2+uj?(k)  ] 
-2ir/a]L  3  ~ 


(2.20) 


To  determine  the  long  wavelength  limit,  we  expand  the  relevant 
expressions  in  terms  of  the  lattice  parameter  aQ  and  keep  only  those 
terms  of  lowest  order.  In  taking  the  limit  as  a  -♦  0,  the  density 

and  the  elastic  constants  must  remain  finite.  This  implies  that 

3  1 

the  mass  m  is  of  order  a  and  the  force  constants  are  of  order  a  . 

o  o 

The  frequencies  of  the  perfect  crystal  normal  modes  in  the 
long  wavelength  limit  can  be  written: 


2,,  x  2,2  2,2  2,2 

u).  (k)  =  c  k  +  c.  k  +  c.  k 
jv~'  jx  x  jy  y  jz  z 


(2.21) 


where  the  coefficients  c  for  each  model  are  found  by  comparison 

J  ** 

with  Equations  (1.12),  (1.13),  and  (1.14)  in  the  longwavelength 


limit.  The  matrix  of  the  values  cf  for  Model  I  is: 

3°/ 


(Model  I) 


6 

P 

3 


3  3 

a  a 

a  a 


(2.22) 


For  Model  II,  we  have  the  added  complication  that  the  frequency 
expressions  for  modes  j  =  2,3  have  a  term  that  involves  the  product 
of  ky  and  k^.  In  the  long  wavelength  limit,  the  normal  mode  fre¬ 


quencies  are: 
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BWJgCk)  =  Xkxao  +  4  (y+X)  kyai  +  4-(y+x>kzai 

~  y X)  ky&^k^a^  (2.23; ) 

mWgOk)"  Xk2a2  +  ^(y+X)k2a2  +  |(y+X)k^ 

+  |(v-X)kya;Lkza1  •  (2.23b) 

These  frequencies  can  be  put  into  the  form  (2.21)  by  a  change  of 
variable. 

J  “2:  kiai  ■  Vi  -  Tmf  \ai 

j  -  3;  kjaj  -  k^j  +  £*=$■  kyaj  (2.24) 

with  the  result: 

«2(y  -  ^yo  +  kM  +  *0+0  0=^) 2 

nuug(k)  =  Xk2a2  +  k2a2  +  l(Y+x)(k^ai)2  (2.25) 

Thus,  if  we  remember  for  the  time  being  that  the  wave  vector  com- 

o 

ponent  k  is  different  for  each  mode  i,  then  the  c.  matrix  for 
z  JQ 

Model  II  is: 


(Model  II) 


(2.26) 
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In  changing  the  integration  variable  from  k  to  k'  or  k"  as 
in  Equations  (2.24),  the  only  effect  is  to  add  a  constant  to  the 
integration  limits.  Before  taking  the  long  wavelength  limit,  the 
integrand  is  periodic  in  the  extended  zone  scheme.  Therefore  add¬ 
ing  a  constant  to  the  limits  of  the  integral  has  no  effect  on  the 
value  of  the  integral.  The  same  is  true  in  the  long  wavelength 
limit  provided  the  new  integration  limits  always  span  the  origin. 
This  last  condition  is  always  met  for  both  cases;  therefore: 

2rr/a,  2rr/a,  2n/a, 


/  dkz  -  f  dk-  -  f  d^ 

-2rt/a,  -2rr/a,  -2rT/a, 


(2.27) 


and  we  can  forget  the  distinction  between  k  ,  k'  ,  and  k"  that 

z  z  z 

o 

we  made  when  we  wrote  the  c .  matrix  in  (2.26) . 

J  ® 
o 

The  J#i/(k  k  i,-y  )  matrix  elements  can  be  written  in  the  long 
v'k  x  y 


wavelength  limit  as: 


Ju(kxkyj,-y  ) 


C,(l) 

- wi - 

4mL 


2rr/a 


dk 


3  (yyM 


z  [d2  +  c2  k2]2 
-2rr/ax  <3  Jz  z 


2  Ci(2) 

J22(kxkyJ ’"y  >  =  “  ^ 


4mL. 


2rr/a, 

/  dk= 
-2n/a, 


<Wa»i 


J12(kxV’'y  5 


[C  .  (1)^(2)]^ 


4mL 


2rr/a, 

/ 

-2TT/a, 


(kz  -  kY)a? 


J*  (k  kj,-y2) 


z  [d2  +  c2  k2]2  21  *  * 

1  3  jz  zJ 


(2.28) 


where : 


,2  2  2,2  2,2 

d  =  y  +  c  k  +  c  k 
J  J  jx  x  jy  y 


(2.29) 


These  expressions  are  valid  for  both  models  and  for  all  j„  The 

limits  of  integration  can  be  extended  to  ±  ®  since  the  integrands 

ai ;  properly  convergent.  The  terms  in  the  integrand  that  are  odd 

in  k  vanish,  and  we  obtain: 
z  1 


Jll(kxkyj,-y2)  =  -C.(1)IV  +  Qky] 


J22(kxV’"y2)  =  -cj(2)£p  +  Qky^l 


1  -iikvai 


J12(kxV'-y  }  =  "[Cj(1)  Cj(2)] 2[P-Qky]  e 


where : 


"  J21<kxV--y2) 


(2.30a) 


P  = 


rra 

4raL 


V-  (cjJ 


Q  - 


_ /_i_\ 

“zVjzU;  / 


(2.30b) 


All  matrix  elements  are  of  order  a^  . 

o 


The  M. . /(k  k  j,-y^)  matrix  elements  in  the  long  wavelength 
x  y 

limit  can  be  written: 


M. ,  (k  k  j,-y  )  = 
11  x  yJ  *  J 


*/2  a  c,(l) 


o  .1 

16nm 


2rr/a 


/ 


-2rr/a. 


1  /,  ,  N  2  2 
(k  +k  )  a, 

dk  -  ~ — — 

z  d'ft-c^  k 
J  J2  z 
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„  J%  a  C  (2)  r^/a  1  (k-k  )2a? 

M22(kxV’'y  ) - T§=i -  /  dk_  -^-z  „?  .  A 


16TTm 


-2n/a . 


z  2  ,  2 

d  ,+c  .  k. 

J  jz  z 


M12(kxV’'y  5  " 


2%  J2  a0[C.(l)C.(2)]s  ijkvaj 
5  IS™  e 


2n/aj  (k2-k2)a2  t 

X  /  dkz  ~2  ~2~  ’  M21(kxkvj’~y  ' 

J  “KV2 

-2TT/a, 


(2.31) 


Again,  these  expressions  are  valid  for  both  models  and  all  j .  Some 
these  integrals  do  not  converge  when  the  limits  of  integration 
are  extended  to  ±  ®.  Hence,  we  retain  finite  limits  for  those 
integrals  that  diverge  and  change  the  limits  of  the  other  integrals 
to  ±  °°  to  obtain: 


Mn<kxV--y> 


CjU)  [R  -  S  +  Tkz] 


M22<kxkyj>'y  ) 


M12(kxkyJ  > 'y2) 


Cj (2)  [R  -  S  +  Tk p 


[C,(1)C.(2)]5  [R  -  S 

u  U 


2 

Tkp  e  y 


Mn1 (k  k  j ,-y  ) 
21  x  yJ *  J  ' 


(2.32) 


where : 


72  agd, 
8mcjz 
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&  a; 


T  = - - — 

8mdjcjz 


(2.33) 


The  coefficient  R  is  of  order  a  *  whereas  S  and  T  are  of  order  a®. 

o  o 

All  three  terras  must  be  retained  in  the  limit  of  a  -*  0,  however, 

o 

2  “1 

until  after  the  inverse  matrix  [I  -  M(kxk^j,-y  ) ] ^ /  is  found. 

o 

The  determinant  A(kxkyj,-y  )  of  the  matrix 
I  -  M(kxkyj,-y2)  is: 

A(kxkyj,-y2)  =  1  -  [Cj(l)+CJ(2))R  +  4^ (DC^ (2)k2 


+  [Cj(l)+Cj(2)][S-Tkp 


-  4TSC,(l)C.(2)k~ 
%}  j  y 


(2.34) 


Of  the  five  terras  on  the  right  hand  side,  the  first  two  are  of 


order  1,  the  second  two  are  of  order  a  ,  and  the  last  is  of  order 
’  o 

2 

a  .  The  two  terms  of  order  1  when  writcen  in  full  are: 
o 


1  -  [Cj(l)+Cj(2)]  R  =  1  -  — 


a„  [C.(l)+C.(2)] 


(2.35) 


These  terms  identically  vanish  because: 

2  4  [C  UHC  (2)] 


(2.36) 


for  all  j  and  for  both  models.  Thus,  the  lowest  order  of  the 

determinant  is  a^  and  can  be  written: 

o 


2  a/2 

Mk  k  j,-y2)  =  - 

y  4c.  d. 

jz  J 


2  2,2  ~2  ,  2  i 

y  +  c  k  +  c  k 
9  JX  X  jy  y 


(2.37) 
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where : 
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o'15  (y2) 


-rr/a. 


-iT/a , 


2  2  ,  2  ~2  ,  2 
y+c  k+c  k 
jx  x  jy  y 


(2.43) 


o  ,2  ~2  .  “1 

nj2)(y2>  “  °  ^  /  ' 


y  d^2 


(2.44) 


-Tt/a  •  “rr/a-] 


(1)  2 

For  the  term  Q:  (y  ) ,  the  limits  for  the  integral  over  k  can  be 

«/ 

extended  to  ±  »  and  the  integration  can  be  carried  out  directly 
to  give: 


na)(y2),_fo_ 

16n8»  L 


x  [y2+c2  k2]^ 
jx  xJ 


(2.45) 


For  this  integral,  the  limits  must  remain  finite  until  after  the 


integration  is  performed. 


fij1>(y2) - rr - }  loe  kAkl+  4 

J  l6m\jycjx  (  I  '  7 


In  the  limit  as  aQ  ■*  0,  this  becomes: 


V^o 


?)•] 


V-n/ao 


0<l)(y2)  -  - 

J 


8ttc  .  c  . 
jx  jy 


log|y|  +  constant 


(2)  2 

For  the  term  Q.  (y  ) ,  we  split  the  integrand  as: 
J 


(2.46) 


(2.47) 


(c^  -c2  ) k2 
d?  d2 


2  2  , 2  ~2  ,  2  2  2,2  2,2 
y  +c  .  k  +c  .  k  y  +c  k  +c  k 
jx  x  jy  y  9  jx  x  jy  y 


(2.48) 
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Thus,  aside  from  different  'coefficients,  the  resulting  two  inte- 
/  o')  2 

grals  for  fi .  (y  )  are  of  exactly  the  same  form  as  the  integral 
( 1) 

for  Q.  (y  )  in  (2.43).  Thus,  we  obtain  directly: 

J 


(y2)  = 


8tt  1  c  .  c  .  c  .  c  . 

\  jx  jy  jx  jy 


logly|  +  constant 


(2.49) 


By  substituting  the  results  (2.47)  and  (2.49)  into  Equation  (2.17) 
and  comparing  with  Equation  (2.2),  we  obtain  for  the  term  A: 


fe  s  — L— 

o3xcJy 


c,  "C. 


(2.50) 


For  Model  I,  the  second  term  in  the  sum  vanishes  because  the 

term  c.  =  c.  for  all  j.  Thus,  by  using  the  identifications  in 

jy  jy  ’ 

(2.21)  and  (2.22)  we  obtain: 


(Model  I) 


8rrrc' 


2  M 

(l+r2)2 


(2.51) 


where : 


(2.52) 


and  c^  and  ct  are  the  longitudinal  and  transverse  velocities  of 
sound  respectively  in  the  bulk  crystal  (100)  directions. 

For  Model  II,  the  second  term  in  (2.50)  vanishes  for  the  mode 
j  ~  1,  and  is  the  same  for  modes  j  =  2  and  j  =  3.  For  this  model, 


the  term  c .  is : 

jy 


j-2,3;  Ejy  -  3  Sj,  - 


(2.53) 


27 


so  that  the  coefficient  A  becomes: 


(Model  II) 


A  = 


!  +  2j±±ihl  ilm  \ 
J2r  l  X  ) 


(2.54) 


g 

Dobrzynski  and  Leman  calculated  the  specific  heat  using 
Model  II  for  the  (100)  surface  by  a  different  method  and  obtained 
the  result: 

(Model  II,  (100)  surface) 

A  -  — ~k  !l  +  2r  |  (2.55) 

8rrr  c^  (  J 


In  the  isotropic  limit  where  r  -*  1  and  x  1,  the  three  results 
(2.51),  (2.54),  and  (2.55)  reduce  tc: 


Aiso 


3S  . 

_ l 

8™l 


(2. bb) 


which  agrees  with  the  result  of  Maradudin  ei. al.1 

The  value  of  the  coefficient  A  for  the  surface  of  an  isotropic 
solid  has  been  determined  from  a  correct  continuum  calculation 
to  be: ^ 

g 

Aiso  =  — ^2  (if)  (continuum)  (2.57) 


The  fact  that  our  value  of  A±so  is  10%  lower  than  the  continuum 
value  is  not  too  surprising  since  our  models  are  not  rotationally 
invariant.  * 


For  small  deviations  from  isotropy,  we  introduce  the  quantity 
A: 


A  =  1  -  r 


(2.58) 


To  lowest  order  in  A,  the  coefficient  A  for  the  three  results 
(2.51),  (2.54),  ,ind  (2.55)  becomes: 

A  =  Aigo  (1  +  -|  A)  (2.59) 

The  fact  that  we  find  the  deviation  of  the  specific  heat  for  small 
deviations  from  isotropy  to  be  independent  of  the  surface  and  model 
for  the  two  surfaces  and  two  models  considered  supports  the  belief 
that  the  expressions  give  meaningful  qualitative  behavior  in  this 
limit. 


III.  CONCLUSIONS 

(1)  We  have  found  the  change  in  the  specific  heat  at  low 
temperatures  due  to  the  (110)  surface  in  a  cubic  monatomic  lattice 
for  two  simple  interatomic  force  constant  models  to  be: 


(Model  I) 

acv(t)  - 


3k|g(3) 

4nft2rc2 


(Model  II) 
ACv(T)  = 


3k|c(3) 

4TTft2rc2 


!  i  +  2 (l±rM 
(  J*  r 


(3.1) 


(3.2) 


To  our  knowledge,  this  is  the  first  determination  of  the  change 
in  the  specific  heat  due  to  a  surface  that  is  anisotropic . 

(2)  We  have  shown  that  the  (110)  surface  of  a  simple  cubic 
monatomic  lattice  gives  rise  to  Rayleigh  waves  that  are  split  off 
from  the  bulk  modes  for  very  simple  force  constant  models.  The 
surface  waves  arise  ever,  for  the  case  where  the  atomic  motion 
normal  to  the  surface  is  totally  uncoupled  from  motion  parallel 
to  the  surface.  The  only  condition  needed  to  obtain  Rayleigh  waves 
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is  that  the  "bonds"  between  atoms  on  opposite  sides  of  the  surface 
be  oblique  to  the  surface, 

(3)  We  have  provided  a  second  example  of  the  practicality 
of  the  Green’s  function  method  suggested  by  Maradulin  and  Ashkin'1® 
in  calculating  the  low  temperature  surface  specific  heat.  The  use¬ 
fulness  of  the  method  depends  upon  the  ability  to  express  the 
resulting  integral  equation  in  separable  form.  Maradudin  and 
Wallis*1  succeeded  in  doing  this  by  assuming  central  force  inter¬ 
actions.  This  paper  succeeded  in  doing  this  by  using  the  fact  that 
the  eigenvectors  for  the  assumed  models  are  independent  of  wave 
vector. 

(4)  We  have  found  that  the  value  of  the  surface  specific 
heat  for  the  two  models  in  the  isotropic  limit  is  10%  lower  tlian 
the  value  found  by  correct  continuum  theory.  We  believe  that  this 
discrepancy  is  due  solely  to  the  fact  that  the  assumed  models  do 
not  satisfy  the  rotational  invariance  condition. 

(5)  We  have  shown  that  as  the  longitudinal  velocity  of  sound 
increases  relative  to  the  transverse  velocity  of  sound  for  values 
that  nearly  satisfy  the  isotropy  condition,  the  value  of  the  surface 
specific  heat  also  increases.  This  result  is  independent  of  the 
model  and  surface  to  first  order  in  the  deviation  from  isotropy 
for  the  two  models  presented  and  for  the  (110)  and  (100)  surfaces. 
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FIGURE  CAPTIONS 


The  atomic  configuration  and  unit  vector  designation 
for  the  (110)  surface  of  a  simple  cubic  monatomic 
lattice.  The  unit  vector  x  is  perpendicular  to  the 
y-z  plane  shown  and  the  atoms  designated  5  and  6  are 

A 

located  at  ±  xa  . 

o 

The  dispersion  curves  for  the  bulk  and  surface  phonons 

A 

propagating  in  the  y  direction  for  Model  II  where  the 
force  constants  are  related  by  y  ~  2X.  The  bulk  mode 
dispersion  curves  for  Model  I  are  the  same  as  those 
shown  here  if  the  force  constants  6  and  p  are  identi¬ 
fied  with  y  and  X  respectively.  With  this  identifica¬ 
tion,  the  surface  mode  j  =  1  is  identical  for  the  two 

models  but  the  degenerate  modes  j  =  2,3  are  higher  for 

1 

Model  I  by  a  factor  of  (9/8) 2 . 

The  relation  between  the  complex  wave  vector  kz  and  the 
wave  vector  k^  for  Models  I  and  II.  The  attenuation  co¬ 
efficients  for  all  modes  in  Model  I  and  for  the  mode 
j  -  1  for  Model  II  are  independent  of  the  force  con¬ 
stants  and  diverge  as  ky  approaches  the  zone  boundary. 
The  real  part  of  the  wave  vector  kz  vanishes  for  these 
modes.  The  real  and  imaginary  parts  of  k  for  modes 
j  =  2,3  in  Model  II  are  shown  for  the  case  y  -  2X. 
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MODEL  I 
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